if(length(grep("smoke$", getwd())) ) {

source("CohortDeffFunction.r")
} else {

source("C:/Users/hjiang/Desktop/sampleSizeGoogle/smoke/CohortDeffFunction.r")

        }
#source("CohortDeffFunction.r")

library(mgcv)
library(MASS)
library(INLA)
library(lme4)
library(glmmBUGS)
library(nlme)
library(gee)
parameters=list(
   size=500,
   DiffEffSizeGroup = 0.75              ,
   DiffEffSizeIndi = 2                   ,
   DiffEffSizeInter = 1.2                 ,
   GroupMissRate = 0.10                     ,
   IndiMissRate = 0.05                       ,
   SmokeMissRate = 0.05                       ,
   probIndiVar = 0.1                           ,
   probGrVar = 0.4                       ,
   regionsd = 0.5    ,
   mu = 4/6
)
Ngroups = c(10, 20, 50, 100)
simulationTime=500


###########################
###### First Senario ######
###########################

varying = NULL
###### effect size = 0 null model
######work nwith .001, .001 #######
nullmodel = T
Nmethods = c("lmer","pql","geeNaive", "geeRobust","inla")
null500 = allEstCI(Ngroups,parameters, varying=NULL,nullmodel,simulationTime, Nmethods )
save(null500, file="null500.Rdata")
null500CovBia = CoverBias( allEstCI=null500, parameters, varying=NULL, nullmodel=T, simulationTime)
SEgraph(allEstCI=null500, CovBia=null500CovBia, nullmodel=T, graphName="null500.pdf")
CoverBiasgraph(allEstCI=null500, Covera=null500CovBia$Covera, 
               biasE=null500CovBia$biasE, 
               nullmodel, 
               CgraphName="Cnull500.pdf", BgraphName="Bnull500.pdf")

####### when low prevalence #######
nullmodel = T
varying = list("mu"=1/9)
Nmethods = c("lmer","pql","geeNaive", "geeRobust","inla")
null500lowpre = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(null500lowpre, file="null500lowpre.Rdata")
null500lowpreCovBia = CoverBias( allEstCI=null500lowpre, parameters, varying=list("mu"=1/9), nullmodel=T, simulationTime)
SEgraph(allEstCI=null500lowpre, CovBia=null500lowpreCovBia, nullmodel=T, graphName="null500lowpre.pdf")
CoverBiasgraph(allEstCI=null500lowpre, Covera=null500lowpreCovBia$Covera, 
               biasE=null500lowpreCovBia$biasE, 
               nullmodel, 
               CgraphName="Cnull500lowpre.pdf", BgraphName="Bnull500lowpre.pdf")

####### when low prevalence and low correlation#######
nullmodel = T
varying = list("mu"=1/9, "regionsd"=0.1)
Nmethods = c("lmer","pql","geeNaive", "geeRobust","inla")
null500lowprelowcor = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(null500lowprelowcor, file="null500lowprelowcor.Rdata")
null500lowprelowcorCovBia = CoverBias( allEstCI=null500lowprelowcor, parameters, varying=list("mu"=1/9,"regionsd"=0.2),
               nullmodel=T, simulationTime)
SEgraph(allEstCI=null500lowprelowcor, CovBia=null500lowprelowcorCovBia, nullmodel=T, graphName="null500lowprelowcor.pdf")
CoverBiasgraph(allEstCI=null500lowprelowcor, Covera=null500lowprelowcorCovBia$Covera, 
               biasE=null500lowprelowcorCovBia$biasE, 
               nullmodel, 
               CgraphName="Cnull500lowprelowcor.pdf", BgraphName="Bnull500lowprelowcor.pdf")


###### when size=5000 ###########
nullmodel= T
varying = list("size"=5000)
Nmethods = c("lmer","pql","geeNaive", "geeRobust","inla")
null5000 = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(null5000, file="null5000.Rdata")
null5000CovBia = CoverBias( allEstCI=null5000, parameters, varying=list("size"=5000), nullmodel=T, simulationTime)
SEgraph(allEstCI=null5000, CovBia=null5000CovBia, nullmodel=T, graphName="null5000.pdf")
CoverBiasgraph(allEstCI=null5000, Covera=null5000CovBia$Covera, 
               biasE=null5000CovBia$biasE, 
               nullmodel, 
               CgraphName="Cnull5000.pdf", BgraphName="Bnull5000.pdf")

####### low prevelence ###########
nullmodel = T
varying = list("size"=5000, "mu"=1/9)
null5000lowpre = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(null5000lowpre, file="null5000lowpre.Rdata")
null5000lowpreCovBia = CoverBias( allEstCI=null5000lowpre, parameters, varying=list("size"=5000, "mu"=1/9), 
                       nullmodel=T, simulationTime)
SEgraph(allEstCI=null5000lowpre, CovBia=null5000lowpreCovBia, nullmodel=T, graphName="null5000lowpre.pdf")
CoverBiasgraph(allEstCI=null5000lowpre, Covera=null5000lowpreCovBia$Covera, 
               biasE=null5000lowpreCovBia$biasE, 
               nullmodel, 
               CgraphName="Cnull5000lowpre.pdf", BgraphName="Bnull5000lowpre.pdf")


####### when low prevalence and low correlation#######
nullmodel = T
varying = list("size"=5000, "mu"=1/9, "regionsd"=0.2)
null5000lowprelowcor = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(null5000lowprelowcor, file="null5000lowprelowcor.Rdata")
null5000lowprelowcorCovBia = CoverBias( allEstCI=null5000lowprelowcor, parameters, 
               varying=list("size"=5000, "mu"=1/9,"regionsd"=0.2), nullmodel=T, simulationTime)
SEgraph(allEstCI=null5000lowprelowcor, CovBia=null5000lowprelowcorCovBia, nullmodel=T, graphName="null5000lowprelowcor.pdf")
CoverBiasgraph(allEstCI=null5000lowprelowcor, Covera=null5000lowprelowcorCovBia$Covera, 
               biasE=null5000lowprelowcorCovBia$biasE, 
               nullmodel, 
               CgraphName="Cnull5000lowprelowcor.pdf", BgraphName="Bnull5000lowprelowcor.pdf")





##############################
###### Second Senario #########
###############################
###### effect size = 0 null model #######
nullmodel = F
varying=NULL
s500 = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime,Nmethods)

save(s500, file="s500.Rdata")
s500CovBia = CoverBias( allEstCI=s500, parameters, varying=NULL, nullmodel=F, simulationTime)
SEgraph(allEstCI=s500, CovBia=s500CovBia, nullmodel=F, graphName="s500.pdf")
CoverBiasgraph(allEstCI=s500, Covera=s500CovBia$Covera, 
               biasE=s500CovBia$biasE, 
               nullmodel, 
               CgraphName="Cs500.pdf", BgraphName="Bs500.pdf")



####### when low prevalence #######
nullmodel = F
varying = list("mu"=1/9)
s500lowpre = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime)
save(s500lowpre, file="s500lowpre.Rdata")
s500lowpreCovBia = CoverBias( allEstCI=s500lowpre, parameters, varying=list("mu"=1/9), nullmodel=F, simulationTime)
SEgraph(allEstCI=s500lowpre, CovBia=s500lowpreCovBia, nullmodel=F, graphName="s500lowpre.pdf")
CoverBiasgraph(allEstCI=s500lowpre, Covera=s500lowpreCovBia$Covera, 
               biasE=s500lowpreCovBia$biasE, 
               nullmodel, 
               CgraphName="Cs500lowpre.pdf", BgraphName="Bs500lowpre.pdf")


####### when low prevalence and low correlation#######
nullmodel = F
varying = list("mu"=1/9, "regionsd"=0.2)

s500lowprelowcor = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(s500lowprelowcor, file="s500lowprelowcor.Rdata")
s500lowprelowcorCovBia = CoverBias( allEstCI=s500lowprelowcor, parameters, varying=list("mu"=1/9,"regionsd"=0.2), 
               nullmodel=F, simulationTime)
SEgraph(allEstCI=s500lowprelowcor, CovBia=s500lowprelowcorCovBia, nullmodel=F, graphName="s500lowprelowcor.pdf")
CoverBiasgraph(allEstCI=s500lowprelowcor, Covera=s500lowprelowcorCovBia$Covera, 
               biasE=s500lowprelowcorCovBia$biasE, 
               nullmodel, 
               CgraphName="Cs500lowprelowcor.pdf", BgraphName="Bs500lowprelowcor.pdf")




###### when size=5000 ###########
nullmodel= F
varying = list("size"=5000)
s5000 = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(s5000, file="s5000.Rdata")
s5000CovBia = CoverBias( allEstCI=s5000, parameters, varying=list("size"=5000), nullmodel=F, simulationTime)
SEgraph(allEstCI=s5000, CovBia=s5000CovBia, nullmodel=F, graphName="s5000.pdf")
CoverBiasgraph(allEstCI=s5000, Covera=s5000CovBia$Covera, 
               biasE=s5000CovBia$biasE, 
               nullmodel, 
               CgraphName="Cs5000.pdf", BgraphName="Bs5000.pdf")


nullmodel = F
varying = list("size"=5000, "mu"=1/9)
s5000lowpre = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(s5000lowpre, file="s5000lowpre.Rdata")
s5000lowpreCovBia = CoverBias( allEstCI=s5000lowpre, parameters, varying=list("size"=5000, "mu"=1/9), nullmodel=F, simulationTime)
SEgraph(allEstCI=s5000lowpre, CovBia=s5000lowpreCovBia, nullmodel=F, graphName="s5000lowpre.pdf")
CoverBiasgraph(allEstCI=s5000lowpre, Covera=s5000lowpreCovBia$Covera, 
               biasE=s5000lowpreCovBia$biasE, 
               nullmodel, 
               CgraphName="Cs5000lowpre.pdf", BgraphName="Bs5000lowpre.pdf")


####### when low prevalence and low correlation#######
nullmodel = F
varying = list("size"=5000, "mu"=1/9, "regionsd"=0.2)
s5000lowprelowcor = allEstCI(Ngroups,parameters, varying,nullmodel,simulationTime, Nmethods)
save(s5000lowprelowcor, file="s5000lowprelowcor.Rdata")
s5000lowprelowcorCovBia = CoverBias( allEstCI=s5000lowprelowcor, parameters, varying=list("size"=5000, "mu"=1/9,"regionsd"=0.2), 
               nullmodel=F, simulationTime)
SEgraph(allEstCI=s5000lowprelowcor, CovBia=s5000lowprelowcorCovBia, nullmodel=F, graphName="s5000lowprelowcor.pdf")
CoverBiasgraph(allEstCI=s5000lowprelowcor, Covera=s5000lowprelowcorCovBia$Covera, 
               biasE=s5000lowprelowcorCovBia$biasE, 
               nullmodel, 
               CgraphName="Cs5000lowprelowcor.pdf", BgraphName="Bs5000lowprelowcor.pdf")
